function output = calculateTfcData_Core(cfg)

[Y,I] = sort(cfg.Input.ListOfTfcFiles);
firstloop = true;
tStart = tic;
for i = I % load in ASCII order...
    disp(['Loading file ' num2str(find(I==i)) '/' num2str(length(I)) '...']);
    load(cfg.Input.ListOfTfcFiles{i});
    disp(['  --> ' data.channel ' of ' data.epochs.name]);
    if firstloop;
        tfc.subject = cfg.SubjectName;
        tfc.epochs = data.epochs.name;
        tfc.times = data.times;
        tfc.freqbins = cfg.frequencies;
        tfc.freqbindescriptions = cfg.freqbindescriptions;
        firstloop = false;
    end
    disp('  --> Calculating...');
    safeChannelName = data.channel;
    safeChannelName(safeChannelName == '-') = '_';
    tfc.data.(safeChannelName) = CreateChannelCollapsedDataByTrialMean(data, cfg.frequencies, cfg);
    disp(['  --> Elapsed time since starting with this epochset: ' num2str(toc(tStart)) 's.']);
end

disp('All channels have been processed.');
disp('  --> Saving...');
filename = [cfg.Output.TfcDataFolder cfg.SubjectPrimaryPrefix ' ' data.epochs.name ' Tfc.mat'];
save(filename, '-v7.3', 'tfc');
output.TfcDatafile = filename;
end

function [tfcollapsed] = CreateChannelCollapsedDataByTrialMean(data, frequencies, cfg)
dabs = abs(data.tf); %TODO waarom hier de absolute waarde nemen? Zou niet nodig moeten zijn?
%% Baseline
switch cfg.BaselineSetup
    case '2sLeft-2sRight'
        % Wel baseline, van -4s tot -2s en van 2s tot 4s...
        times = data.times;
        baselineTf = (times>-4000 & times<-2000) | (times>2000 & times<4000);
        baselinePerTrial = mean(dabs(:,baselineTf,:),2);
        dabsm = dabs ./ repmat( baselinePerTrial, [1,size(dabs,2),1]);
    case 'New1s'
        stimulus = ~isempty(strfind(data.epochs.name, 'Stimulus'));
        includeEpochs = true(size(data.markers));
        
        if stimulus
            times = data.times;
            baselineTf = (times>-1000 & times<0);
            baselinePerTrial = mean(dabs(:,baselineTf,:),2);
            if cfg.RemoveEpochsConditionally
                for epI = 1:length(data.markers)
                    latMs = [data.markers{epI}.latencyMs];
%                     includeEpochs(epI) = isempty(find(latMs>-3000 & latMs<0));

                    % No events allowed
                    area1Ok = isempty(find(latMs>-2000 & latMs<0, 1));
                    % Only one stim. event allowed
                    sss = find(latMs>0 & latMs<2000);
                    area2Ok = isempty(sss) || ((length(sss) == 1) && (data.markers{epI}(sss).type(1) == 'R'));
                    % Update include-matrix
                    includeEpochs(epI) = area1Ok & area2Ok;
                end
            end
            
        else
            times = data.times;
            baselineTf = (times>-3000 & times<-2000);
            baselinePerTrial = mean(dabs(:,baselineTf,:),2);
            if cfg.RemoveEpochsConditionally
                for epI = 1:length(data.markers)
                    latMs = [data.markers{epI}.latencyMs];
                    % No events allowed
                    area1Ok = isempty(find(latMs>-3000 & latMs<-2000, 1));
                    % No events allowed
                    area2Ok = isempty(find(latMs>0 & latMs<1000, 1));
                    % Only one stim. event allowed
                    sss = find(latMs>=-2000 & latMs<0, 1);
                    area3Ok = isempty(sss) || ((length(sss) == 1) && (data.markers{epI}(sss).type(1) == 'S'));
                    % Update include-matrix
                    includeEpochs(epI) = area1Ok & area2Ok & area3Ok; % 
                    area(1,epI) = area1Ok;
                    area(2,epI) = area2Ok;
                    area(3,epI) = area3Ok;
                end
                disp(area);
            end
        end

        dabsm = dabs;% ./ repmat( baselinePerTrial, [1,size(dabs,2),1]);
        if cfg.RemoveEpochsConditionally
            dabsm = dabsm(:,:,includeEpochs);
            disp(['  --> Removing ' num2str(sum(~includeEpochs)) ' of ' num2str(length(includeEpochs)) ' epoch(s)...']);
        end
    otherwise
        % Geen baseline...
        dabsm = dabs ./ repmat( mean(dabs, 2), [1,size(dabs,2),1]); %Dit is geen normalisatie van variantie, maar alleen van gemiddelde. Overleg met Maartje.
end

%% Collapsing
for i=1:length(frequencies)
    tfcollapsed(i).frequencyrange = frequencies(i,:);
    bin = dabsm(data.freqs>=frequencies(i,1) & data.freqs<=frequencies(i,2),:,:); % Gets 'freq in bin' x 'time' x 'trials'
    timesByTrials = squeeze(mean(bin,1)); %calculate average per time point per trial for each freq bin
    tfcollapsed(i).mean = mean(timesByTrials,2); % calculate average over all trials for each time point
    tfcollapsed(i).std = std(timesByTrials,[],2);
    tfcollapsed(i).qEpochs = size(dabsm,3);
end

end